//options
clear all
set maxvar  32000, permanently
set matsize 11000, permanently
set more off, permanently
global bootstraps 1000
set seed 23

// macros
global klmshare				: env klmshare
global projects				: env projects
global klmperry             : env klmperry
global storageb				: env storageb

global perrydatas       = "${klmperry}/CBA/data/perry/raw"
global perrydata        = "${klmperry}/CBA/data/perry/clean"
global masterinter  	= "${klmshare}/CurrentRAs/FredB/test"
global output           = "${projects}/ece_parenting/tex_files/other"
global dataperry        = "$klmperry/PerryPreschool/Data/Perry_PARI_and_Other_Data/PARI/DATA_PARI/"
global abcjpeanalysis   = "${klmshare}/Data_Central/Abecedarian/data/ABC-CARE/extensions/cba-iv/"
global dataihdp         = "${projects}/ece_parenting/data"
global dataihdp         = "$storageb/dc_data/"

cd  $dataihdp
use ihdp_data, clear
rename _all, lower

replace   bw = bw/1000

// impute iq
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.
foreach var of varlist stndscor iqcage { 
	replace `var' = mdi_24cor    if stndscor ==. & iqcage ==.
}
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.

// low birthweight
gen lb = 1
rename tot_siblings_natural siblings_natural

// list relevant variables
global baseline_child      twin sex bw lb anga black hispanic 
global baseline_mother     mage meduc works married 
global baseline_household  welfare siblings_natural employed_adult
global baseline_economy    employment medinc gpc
global outputs             iqcage stndscor 
global inputs              cum_avg_daycare_36m_sum

// mark baseline sample
reg  $baseline_child $baseline_mother $baseline_household $baseline_economy $outputs $inputs
gen  sample1 = e(sample) 
keep if sample1 == 1

// parenting latents, ages 1 and 3
foreach var of varlist alt_subscale_1_12-alt_subscale_6_12 alt_subscale_1_36-alt_subscale_8_36 {
	summ `var'
	gen  `var'_std = (`var' - r(mean))/r(sd)
}
# delimit
sem (_cons@0 X -> alt_subscale_1_12_std) (_cons@0 X -> alt_subscale_2_12_std) (_cons@0 X -> alt_subscale_3_12_std) 
	(_cons@0 X -> alt_subscale_4_12_std) (_cons@0 X -> alt_subscale_5_12_std) (_cons@0 X -> alt_subscale_6_12_std), method(adf);
predict parenting_age1, latent;
sem (_cons@0 X -> alt_subscale_1_36_std) (_cons@0 X -> alt_subscale_2_36_std) (_cons@0 X -> alt_subscale_3_36_std) (_cons@0 X -> alt_subscale_4_36_std)
	(_cons@0 X -> alt_subscale_5_36_std) (_cons@0 X -> alt_subscale_6_36_std) (_cons@0 X -> alt_subscale_7_36_std) (_cons@0 X -> alt_subscale_8_36_std), 
																																						 cov(e.alt_subscale_1_36_std*e.alt_subscale_2_36_std)
																																						 cov(e.alt_subscale_3_36_std*e.alt_subscale_4_36_std)
																																						 cov(e.alt_subscale_5_36_std*e.alt_subscale_6_36_std)
																																						 cov(e.alt_subscale_7_36_std*e.alt_subscale_8_36_std)
																																						 cov(e.alt_subscale_4_36_std*e.alt_subscale_7_36_std)
																																						 cov(e.alt_subscale_2_36_std*e.alt_subscale_7_36_std) method(adf);		
 predict parenting_age3, latent;
# delimit cr
egen    parenting_ages13 = rowmean(parenting_age1 parenting_age3)

// standardize latents in sample
foreach var of varlist parenting_ages13 {
	summ    `var' if tg == 0
	gen     `var'_std = (`var' - r(mean))/r(sd)
	replace `var'_std = `var'_std + 100
}
// log and standardize childcare variable 
summ    cum_avg_daycare_36m_sum  if tg == 0 
replace cum_avg_daycare_36m_sum  = (cum_avg_daycare_36m_sum    - r(mean))/r(sd)
replace cum_avg_daycare_36m_sum  =  cum_avg_daycare_36m_sum    + r(mean)

rename cum_avg_daycare_36m_sum c
rename parenting_ages13_std    p

global inputs_std c p
keep   ihdp site tg bwg $baseline_child $baseline_mother $baseline_household $baseline_economy $inputs_std $outputs

// describe by sample
global sample_all          if twin != .
global sample_singletons   if twin == 0
global sample_twins        if twin == 1

cd $output
foreach group in all singletons twins {
	// treatment effects
	foreach var of varlist $inputs_std $outputs {
		bootstrap, reps($bootstraps) strata(site tg bwg) saving(tebd_`var'_`group', replace): reg `var' tg ${sample_`group'}
		matrix te_`var'_`group' = e(b)[1,1]
		preserve
		use tebd_`var'_`group', clear
		keep _b_tg
		mkmat *, matrix(tebd_`var'_`group')
		matrix colnames tebd_`var'_`group' = tebd_`var'_`group'
		restore
		svmat tebd_`var'_`group', names(col)
		erase tebd_`var'_`group'.dta
	}
	// production functions
	foreach var of varlist $outputs {
		bootstrap, reps($bootstraps) strata(site tg bwg) saving(cpbd_`var'_`group', replace): reg `var' c p
		matrix c_`var'_`group' = e(b)[1,1]
		matrix p_`var'_`group' = e(b)[1,2]
		preserve
		use cpbd_`var'_`group', clear
		keep _b_c _b_p
		rename _b_c cbd_`var'_`group'
		rename _b_p pbd_`var'_`group'
		mkmat *, matrix(cpbd_`var'_`group')
		matrix colnames cpbd_`var'_`group' = cbd_`var'_`group' pbd_`var'_`group'
		restore
		svmat cpbd_`var'_`group', names(col) 
		erase cpbd_`var'_`group'.dta
	}
	// decomposition
	foreach var of varlist $outputs {
		foreach vary in $inputs_std {
			matrix `vary'd_`var'_`group'   = `vary'_`var'_`group'*te_`vary'_`group'
			gen    `vary'dbd_`var'_`group' = `vary'bd_`var'_`group'*tebd_`vary'_`group'
		}
		matrix resd_`var'_`group'   = te_`var'_`group'   - (cd_`var'_`group' + pd_`var'_`group') 
		gen    resdbd_`var'_`group' = tebd_`var'_`group' - (cdbd_`var'_`group' + pdbd_`var'_`group')
	}
}

// keep variable
keep tebd_c_all-resdbd_stndscor_twins
keep if tebd_c_all !=.
collapse (sd) *

foreach group in all singletons twins {
	foreach var in $inputs_std $outputs {
		matrix colnames te_`var'_`group' = te_`var'_`group'
		svmat  te_`var'_`group', names(col)
	}
	foreach var in $outputs {
		foreach vary in $inputs_std res {
			matrix colnames `vary'd_`var'_`group' = `vary'd_`var'_`group'
			svmat  `vary'd_`var'_`group', names(col)
			
		}
	}
}
aorder 

// plot
foreach num of numlist 0(1)20 {
	gen n`num'  =  `num'
	gen nn`num' =  `num' + .5
	gen m`num'  = -`num'
	gen mm`num' = -`num' - .5
}
gen zero = 0


foreach var in all twins singletons {
	foreach vary in iqcage stndscor {
		foreach varyy in cd pd resd {
			gen     `varyy'bdt_`vary'_`var' = abs(`varyy'_`vary'_`var' / `varyy'bd_`vary'_`var')
			gen     `varyy'bdp_`vary'_`var' = (1 - normal(`varyy'bdt_`vary'_`var'))
			drop    `varyy'bd_`vary'_`var' `varyy'bdt_`vary'_`var'
			rename  `varyy'bdp_`vary'_`var' `varyy'bd_`vary'_`var'  
		}
	}
}

foreach var in all twins singletons {
	foreach varyy in c p {
		gen    tebdt_`varyy'_`var' = abs(te_`varyy'_`var' / tebd_`varyy'_`var')
		gen    tebdp_`varyy'_`var' = 1 - normal(tebdt_`varyy'_`var')
		drop   tebd_`varyy'_`var' tebdt_`varyy'_`var'
		rename tebdp_`varyy'_`var' tebd_`varyy'_`var'
	}
}

foreach var in all twins singletons {
	foreach varyy in iqcage stndscor {
		gen    tebdt_`varyy'_`var' = abs(te_`varyy'_`var' / tebd_`varyy'_`var')
		gen    tebdp_`varyy'_`var' = 2*(1 - normal(tebdt_`varyy'_`var'))
		drop   tebd_`varyy'_`var' tebdt_`varyy'_`var'
		rename tebdp_`varyy'_`var' tebd_`varyy'_`var'
	}
}

foreach group in all singletons twins {
	foreach var in $outputs {
		tostring   te_`var'_`group', gen(te_`var'_`group's)   force format(%15.2fc)
		tostring tebd_`var'_`group', gen(tebd_`var'_`group's) force format(%15.2fc)	
		replace te_`var'_`group's = te_`var'_`group's + "({it:p} = " + tebd_`var'_`group's +")"
	
		foreach vary in $inputs_std res {
			tostring   `vary'd_`var'_`group', gen(`vary'd_`var'_`group's)   force format(%15.2fc)
			tostring `vary'dbd_`var'_`group', gen(`vary'dbd_`var'_`group's) force format(%15.2fc)	
			replace `vary'd_`var'_`group's = `vary'd_`var'_`group's + "({it:p} = " + `vary'dbd_`var'_`group's +")"
		}
	}
}

// shift for display
foreach var in $outputs {
	replace pd_`var'_singletons   = pd_`var'_singletons + cd_`var'_singletons
	replace resd_`var'_singletons = pd_`var'_singletons + resd_`var'_singletons
}

foreach var in $outputs {
	replace resd_`var'_twins      = resd_`var'_twins    + cd_`var'_twins
}

cd $output
#delimit
twoway (rbar te_iqcage_singletons zero n9                  , horizontal color(gs0)  yline(5.5, lwidth(vvthick) lpattern(solid) lcolor(gs14)))
	   (rbar cd_iqcage_singletons zero n7                  , horizontal color(gs12) xline(0,   lwidth(vvthick) lpattern(solid) lcolor(gs14)))
	   (rbar cd_iqcage_singletons pd_iqcage_singletons n7  , horizontal color(gs6)  xline(-5,  lpattern(solid) lcolor(gs14)))
	   (rbar resd_iqcage_singletons pd_iqcage_singletons n7, horizontal lcolor(gs0) fcolor(none))
	   
	   (scatter n10 n1                 , msymbol(none) mlabel(te_iqcage_singletonss)   mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter n8  nn1                 , msymbol(none) mlabel(cd_iqcage_singletonss)   mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter n8  nn5                 , msymbol(none) mlabel(pd_iqcage_singletonss)   mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter n8  nn8                 , msymbol(none) mlabel(resd_iqcage_singletonss) mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   
	   
	   (rbar te_iqcage_twins zero n4                , horizontal color(gs0))
	   (rbar cd_iqcage_twins zero n2                , horizontal color(gs12))
	   (rbar zero pd_iqcage_twins                 n2, horizontal color(gs6))
	   (rbar cd_iqcage_twins resd_iqcage_twins    n2, horizontal lcolor(gs0) fcolor(none))
	   
	   (scatter n5 n1                   , msymbol(none) mlabel(te_iqcage_twinss)   mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter n3  nn1                 , msymbol(none) mlabel(cd_iqcage_twinss)   mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   
	   (scatter n3  mm1                 , msymbol(none) mlabel(pd_iqcage_twinss)   mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter n3  nn5                 , msymbol(none) mlabel(resd_iqcage_twinss) mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   
	   ,
		  legend(label(1 "{&Delta}: Treatment Effect") label(2 "Due to Childcare")
											           label(3 "Due to Parenting")
											           label(4 "Residual")
			cols(2) rows(2) order(1 2 3 4) size(medsmall) position(6) region(lcolor(gs0)))
				 
		  ylabel(10 " " 8 "Singletons" 3 "Twins" 1 " ", labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
		  xlabel(-5[5]10,  labsize(medsmall) grid glcolor(gs14) glpattern(solid))
		  ytitle("", size(medsmall))
		  xtitle("Test Score Points")
		  graphregion(color(white)) plotregion(fcolor(white)); 
#delimit cr
graph export iqcage_decomposition.eps, replace

#delimit
twoway (rbar te_stndscor_singletons zero n9                  , horizontal color(gs0)    yline(5.5, lwidth(vvthick) lpattern(solid) lcolor(gs14)))
	   (rbar cd_stndscor_singletons zero n7                  , horizontal color(gs12)   xline(0,   lwidth(vvthick) lpattern(solid) lcolor(gs14)))
	   (rbar cd_stndscor_singletons pd_stndscor_singletons n7  , horizontal color(gs6)  xline(-5,  lpattern(solid) lcolor(gs14)))
	   (rbar pd_stndscor_singletons resd_stndscor_singletons n7, horizontal lcolor(gs0) fcolor(none))
	   
	   (scatter n10 n1                  , msymbol(none) mlabel(te_stndscor_singletonss)   mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter n8  nn1                 , msymbol(none) mlabel(cd_stndscor_singletonss)   mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter n8  n5                 , msymbol(none) mlabel(pd_stndscor_singletonss)   mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter n8  nn7                 , msymbol(none) mlabel(resd_stndscor_singletonss) mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   
	   
	   (rbar te_stndscor_twins zero n4                , horizontal color(gs0))
	   (rbar cd_stndscor_twins zero n2                , horizontal color(gs12))
	   (rbar zero pd_stndscor_twins                 n2, horizontal color(gs6))
	   (rbar cd_stndscor_twins resd_stndscor_twins    n2, horizontal lcolor(gs0) fcolor(none))
	   
	   (scatter n5 n1                   , msymbol(none) mlabel(te_stndscor_twinss)   mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter n3  nn1                 , msymbol(none) mlabel(cd_stndscor_twinss)   mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   
	   (scatter n3  mm1                 , msymbol(none) mlabel(pd_stndscor_twinss)   mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   (scatter n3  n6                 , msymbol(none) mlabel(resd_stndscor_twinss) mlabposition(6) mlabcolor(gs0) mlabsize(vsmall))
	   
	   ,
		  legend(label(1 "{&Delta}: Treatment Effect") label(2 "Due to Childcare")
											           label(3 "Due to Parenting")
											           label(4 "Residual")
			cols(2) rows(2) order(1 2 3 4) size(medsmall) position(6) region(lcolor(gs0)))
				 
		  ylabel(10 " " 8 "Singletons" 3 "Twins" 1 " ", labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
		  xlabel(-5[5]10,  labsize(medsmall) grid glcolor(gs14) glpattern(solid))
		  ytitle("", size(medsmall))
		  xtitle("Test-Score Points")
		  graphregion(color(white)) plotregion(fcolor(white)); 
#delimit cr
graph export stndscor_decomposition.eps, replace
